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o: 

We develop a dynamical theory of heat transfer between two nano systems. In 
| particular, we consider the resonant heat transfer between two nanoparticles due to 

■ the coupling of localized surface modes having a finite spectral width. We model the 

coupled nanosystem by two coupled quantum mechanical oscillators, each interacting 
with its own heat bath, and obtain a master equation for the dynamics of heat 
*^ transfer. The damping rates in the master equation are related to the lifetimes of 

D ■ localized plasmons in the nanoparticles. We study the dynamics towards the steady 



X 



state and establish connection with the standard theory of heat transfer in steady 
state. For strongly coupled nano particles we predict Rabi oscillations in the mean 



■ occupation number of surface plasmons in each nano particle. 
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I. INTRODUCTION 

The theory of heat transfer between nanosystems has been developed for a steady state 
which is usually characterized by time scales much larger than the relaxation time scale of 
surface plasmon excitations. However, such time scales are nowadays accessible experimen- 
tally In particular, in the case of long range plasmons {^J these relaxation times could 
be quite long, demanding the study of the dynamics of heat transfer. 

In addition, nearly all the works studying radiative heat transfer at the nanoscale rely 
on Rytov's flucutational elextrodynamics |3| which is based on the fluctuation-dissipation 
theorem of the second kind jj, 0], and on macroscopic Maxwell equations. This theoretical 
framework was very succesful in the past. Indeed, the seminal works by Lifshitz, and by 
Polder and van Hove , paved the way for numerous studies of Casimir-Lifshitz forces 
and nanoscale radiative heat transfer. As a matter of fact, the main assumption in Rytov's 
theory is that the considered systems are in local thermal equilibrium. By introducing 
macroscopic thermal fluctuating currents or dipole moments it is therefore possible to relate 
the correlation functions of the currents or dipoles to the material properties by applying 
the fluctuation-dissipation theorem. Therefore, the theory is a phenomenological one which 
has no microscopic justification and enables one to study stationary or quasi-stationary 
situations only [8(. A genuinely microscopic quantum mechanical description for the heat 
transfer problem was provided by Janowicz et al. j^] using the Caldera-Legget model. But 
again, this work restricts itself to the steady state, showing that this quantum mechanical 
description leads to the same results for the heat transfer problem as Rytov's theory. 

In this work, we will introduce a simple quantum mechanical model which allows us to 
study the dynamics of the heat transfer rate between two nanosystems in general and be- 
tween two nanoparticles in particular. To keep our model simple we will right from the start 
assume that we have only two nanoparticles supporting localized surface modes which are 
assumed to provide the main heat flux channel [loj |. By this assumption we neglect any 



contribution due to eddy currents 111. Il2fl. crossed electromagnetic interaction 13], multi- 
poles [3], non-Debye relaxation inside the nanoparticle [l5], and many-particle effects 16]. 
Then, the nanoparticles can be modeled as quantum mechanical oscillators which are cou- 
pled due to the interaction of the localized surface modes through their electromagnetic 
fields. Since we are interested in the heat flow, both nanoparticles are coupled to their own 
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heat baths. By deriving the master equation for this system we can determine the dynamics 
of the heat transfer rate towards the steady state. We establish the connection of our model 
to the results obtained from Rytov's theory [ll, 17, 18] by comparison with the steady state 
solutions of our model. 

The structure of the paper is as follows: In section II we motivate and introduce our 
model. The master equation is derived in section III, where we also determine the analytical 
solutions as well as the short- and large-time limits of these solutions and the lowest order 
perturbation results. In section IV we define the heat transfer rate and show how it can be 
derived using Fermi's golden rule. Finally, in section V we compare the steady-state result 
of our model with the known steady-state solutions for the heat transfer rate between two 
nanoparticles determining the coupling constant parameter of our model. In this section we 
also discuss the relaxation dynamics of the heat transfer rate between two nanoparticles. 



II. PHYSICAL MODEL FOR THE DYNAMICS OF HEAT TRANSFER 

BETWEEN NANO SYSTEMS 

The aim of our work is to study the heat transfer dynamics between two nanosystems. 
In particular, we model the resonant heat transfer between two nanoparticles which support 
localized surface modes by means of a quantum mechanical approach. Our model is moti- 
vated by the fact that nanoparticles or nanosystems which have an extension smaller than 
the wavelength of the impinging electromagnetic field can be described by a dipole moment 
p which is induced by the field. This induced dipole moment is related to the incoming field 
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p = e a{u)E (1) 

where a is the polarizability of the nanoparticle and eo is the permittivity of vacuum. The 
resonance of the nanoparticle is determined by the poles of the polarizability. For metallic 
nanoparticles this resonance can be attributed to collective charge density oscillations within 
the particle which are called localized surface plasmons. 

Considering a very small spherical nanoparticle the polarizability can be expressed as [ijl 

where r is the radius of the nanoparticle and e is its permittivity. From this expression 
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for the polarizability it becomes apparent that the resonance of the nanoparticle is given 
by e(u) = —2. This implicit equation can in general be solved for a complex frequency 
oj = u sp — ir which determines the oscillation frequency u sp of the localized charge oscillation 
and its damping or spectral width T due to losses inside the nanoparticle. For Drude metals 
with the permittivity 

e(u) = eoo — ^— (3) 

where u p is the plasma frequency, e OT the background permittivity and 7 a phenomenological 
damping constant we find u sp = u p / V e oo + 2 and r = 7/2 for uj p ^> 7 which is in general 
the case. 
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Figure 1: Sketch of the situation considered. Two nanoparticles separated by a distance d are 
effectively replaced by two dipoles oscillating with the surface mode resonance frequency w sp . 

Let us now assume that we have two identical nanoparticles A and B as sketched in Fig. [TJ 
Particle A has some nonzero temperature 7\ and particle B has temperature T 2 = 0. Then 
localized surface modes will be excited inside particle A due to the thermal fluctuations of 
the charges inside that particle if K-qTi « hu sp , where Kb is Boltzmann's constant and fi is 
Planck's constant. The thermally excited resonant charge oscillations generate a dipole field 
which will excite localized charge oscillations inside particle B. Since the charge oscillations 
are damped inside particle B a part of the exitation energy will be converted into heat and 
particle B heats up, whereas particle A cools down. Of course, the charge oscillations inside 
particle B will also excite charge oscillations in particle A and so forth. So that we can 
expect that energy will be transferred back and forth between the particles until a steady 
state or equilibrium situation is reached. 

The heat flux between two nanoparticles within the steady state can be calculated by 

n n n 

using the fluctuation-dissipation theorem [4j, |5[ and was for example studied in Ref. [11]. 
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Here, we are interested in the dynamics of the relaxation towards the steady state which 
amounts to apply a full non-equilibrium description. To this end, we first chose to describe 
the two nanoparticles as quantum mechanical oscillators A and B having a frequency u sp . 
The Hamiltonian of the two oscillators is 

Hq = hu^a^a + hu S ptfb. (4) 

Here a', a and b\ b are the creation and annihilation operators associated to oscillator A 
and B which fulfill the bosonic commutation relations [a, a* 

] = [6, ftt] = i. Note, that the 
zero point contribution is not taken into account, since it does not affect the dynamics. 

The interaction of the two nanoparticles which is attributed to the coupled charge os- 
cillations is now taken into account by a linear coupling between oscillator A and B of the 
form 

Hi = hg(tfa + a%). (5) 

That means we allow for the exchange of photons (single photon process) between the two 
oscillators. The first term describes the photon transfer from particle A to particle B, i.e. 
|^a,^b) - > |^a — 1,^b + 1), whereas the second term describes the photon transfer from 
particle B to particle A, i.e. \n^,n-B) — > |jia + 1 ; ^b — 1)- The quantity g can be estimated 
by a variety of methods — one of which is described in Sec V. The most direct and general 



method is to use an analog of the mode-mode coupling theory as in case of waveguides [20]. 
The parameter g is essentially given by the overlap of the plasmon field produced by the 
nanoparticle A with the plasmon field of the nano particle B. Note that ii g ^ T then we 
are in the perturbative regime, however for g > T we have non-perturbative regime where 



Rabi oscillations occur. The situation is somewhat reminiscent of cavity QED 21] 



Since both nanoparticles can attain a temperature we couple each oscillator to its own 



251 ] described by the two 



heat bath which consist of a spectrum of independent oscillators 
Hamiltonians 

H B1 = ^2 l huj lj a\ j a lj , (6) 

3 

3 

Again the creation and annihilation operators a\j, a\j and a 2 j, a>2j fulfill the bosonic com- 
mutation relations [aij,ay = [a2i,aL] = 6ij] u\j and w^j are the oscillator frequencies. We 
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Figure 2: Sketch of our model consisting of two harmonic oscillator A and B interacting through Hj. 
Each oscillator is coupled to a heat bath of independend oscillators in equillibrium at temperatures 
T\ and T 2 . 

assume that both heat baths are in thermal equilibrium at temperatures T\ and T 2 , i.e. the 
density operators of the heat baths are given by 

g — /3l/2^Bl/B2 

^ B1 /B2 = ^j e - ft/s /f B1/BJ j (8) 

where (3i/ 2 = 1/KbTi/ 2 is the inverse temperature. 

We further assume that each nanoparticle described by the harmonic oscillators A and B 
is linearly coupled to its heat bath Bl and B2. The corresponding Hamiltonians are given 
by 

Ha-bi = fc/.gijia + a T )(aij - a\j)> ( 9 ) 
3 

H B ^ B2 = faJ2 92 ^ b + ~ °2i) ( 10 ) 

3 

introducing the coupling strengths g\j and g 2 j which are related to the damping T of the 
charge oscillations modeled by the oscillators A and B. When assuming that both reservoirs 
are in thermal equilibrium at temperatures T\ and T 2 and there is no interaction between 
oscillator A and B (g = 0), then by means of this coupling to their heat baths both oscillators 
will in the long-time limit reach an equilibrium state described by the reduced density 
operators 

pA / B = Tr (e-^A/B) 

where = hoj sp a^a and = huj sp b'b. Hence, the two harmonic oscillators A and B would 
acquire the temperatures of their heat baths if there is no coupling, i.e. g = 0. Since in our 
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model the coupling g ^ between the oscillators A and B they will in general not acquire 
the temperatures of their reservoirs in the long-time limit. 

Putting all the different terms together the Hamiltonian of our model sketched in Fig. [2] 
is given by 

H = H + Hi + Hbi + Hb2 + Ha~bi + -^b-B2- (12) 

For later reference we introduce H$ = H + Hi which is the Hamiltonian of a system of two 
coupled harmonic oscillors without any reservoir. 



III. MASTER EQUATION 

In order to determine the dynamics of the system of two coupled harmonic oscillators 
we first derive the master equation of the reduced coupled-oscillator system H$. Using the 
standard Born-Markov approximation in combination with the rotating wave approxima- 



tion 



25[ and tracing out the bath variables using the density operators pei and pB2 from 



Eq. ([8]) we obtain the master equation 

dps . r t 



(13) 



-iu a [a t a,p s ] -iu b [tfb,p s ] 

- \g[a ] b + b j a,p s ] 

— Ki(ni + l)(a)aps — 2aps<x [ + Ps a ^ a ) 

— KirTi (aa) ps — 2a) ps<x + ps^o)) 

- K 2 (n 2 + 1) (tfbps - 2bp s b ] + p s b ] b) 

(btfps - 2tfp s b + psbtf) 

where the coupling constants 

«i = -K^glfi^j -uy), (14) 

3 

K 2 = 7T ^ 9 2 j$(u2j ~ Ub) (15) 
j 

can in our model be identified as the linewidths T of the localized surface modes. We have 
also introduced the mean occupation numbers 

= e n h J a/b _ I (16) 



s 



stemming from the trace over the bath oscillators which are assumed to be in thermal 
equillibrium at temperatures 7\ and T 2 . For the sake of generality and later use we have 
here assumed that the frequencies u a and u b of the oscillators A and B are different. In a 
later stage, we will set these frequencies to u a = u b = w sp . 

The master equation allows us now to determine the dynamical equations of the mean 
values (a^a), (b^b), (a^b) and (tfa). By using the commutation relations of the creation and 
annihilation operators for both oscillators we obtain the set of equations 



_( a t a ^ = -ig^b) - (tfa)) - 2K 1 (a t a) + 2«ini, (17) 

^-(tfb) = -ig{(tfa) - (a+6» - 2n 2 (tfb) + 2n 2 n 2 , (18) 

^(fetfl) = n ab (tfa) - ig({tfb) - (ata)), (19) 

^(atfe) = tt ba (a^b) - \g{^a) - (tfb)), (20) 

where for the sake of clarity we have introduced the new quantities 

Q ab = -i(u a - u b ) - «i - « 2 , (21) 

Q ba = +i(uj a - u b ) - K X - k 2 . (22) 



A. Steady state solutions 

The steady state solutions of the dynamical equations can now be obtained by setting all 
time derivatives ^(a*a), ^{b^b) etc. equal to zero. Solving the resulting set of equations we 
obtain 

+ A^ni + K 2 n 2 ) + 2KiK 2 ni 

(a a) = t, ; , 23 

X 1 A{k x + K 2 ) + 2kxk 2 v ; 

/, fn _ AfjWh + «2W 2 ) + 2/tiK 2 n 2 , s 

1 ' ~ + « 2 ) + 2kik 2 ' 1 j 

(6t a > = -f 2 «i«^-"i) , (25 ) 
(i6 (ki + k 2 )A+2k 1 k 2 ' 

(atfo) = -^a), (26) 



with 



2 ^a6 + ^fea 2 2(^1 + K 2 ) 



ilab^ba (Ki + K 2 ) 2 + (CJ a - CU 6 ) 2 
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The steady state expressions allow us to check the plausibility of our approach. If we set 
for example the coupling between the two oscillators to zero (g = 0), then we find (a^a) = n\ 
and (b^b) = n 2 . That means that the two uncoupled oscillators acquire the temperatures of 
their heat baths, as expected. If we cut off one of the heat baths by setting k\ = (k 2 = 0), 
then we find (a^a) = (tfb) = n 2 ((a 'a) = (tfb) = ni) which means that in this case both 
oscillators take the temperature of the remaining heat bath. Finally, if we let the coupling 
between the two oscillators go to infinity (/ — > oo) then we find 

(a)a) = (tfb) = Kl n x + " 2 n 2 . (28) 

K\ + K 2 «4 + K 2 

In this case, both oscillators are not in an equilibrium state with one of the reservoirs 
anymore. Their mean occupation number is the sum of the equilibrium occupation numbers 
of both reservoirs weighted by the relative coupling strength. 



B. Full dynamical solutions 



We will now derive the full dynamical solutions by setting u a = cj& = cu sp . This corre- 
sponds to the situation of a resonant coupling in which we are interested and it simplifies 
our problem, since then Q a f, = Q^a = —(k>i + ^2)- It follows that the set of four coupled 
differential equations can be recasted into a set of only 3 coupled differential equations by 
substracting the last two differential equations (fT9"|) and (T2"U|) from each other and by con- 
sidering (b^a) — (a^b) as a new dynamical variable. Furthermore we assume that the second 
heat bath has zero temperature, i.e. we set n 2 = 0. Then the 3 coupled differential equations 
can be stated as 

x = Ax + a (29) 
with x = ((ala), (tfb), (tfa) — (atb)) 1 , a = (2ftirTi, 0, 0), and 



/ 



A 



9 



\ 



(30) 



-2k ± 

-2k 2 -g 
\ 2 9 -2g -(ki + k 2 ) J 

This set of coupled first order differential equations can be solved by hand in a standard 
fashion. First we determine the eigenvalues of A by solving the characteristic equation 
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det(A — Al) = 0. We obtain the three eigenvalues 



Ai = + « 2 ), 

A 2 ,3 = -(«i + K 2 ) ± - K2) 2 + V. 



(31) 
(32) 



Then we diagonalize matrix A by introducing the matrix S consisting of the eigenvectors of 
A such that 



^Ai 0^ 



S _1 AS 



V 



B. 



(33) 



A 2 

A 3 y 

By means of the matrices S and S _1 we decouple the set of coupled differential equations 
which translate into 

y = By + b. (34) 

Here y = $ _1 x and b = $ _1 a. Since a = (2K{n\, 0, 0)*, b equals the first column of S _1 
times 2K\fii. We have now three decoupled inhomogeneous differential equations (i = 1, 2, 3) 



which have the solution 



Vi 



y% = \yi + k 



^% + (e Al * - 1) 



X, 



(35) 



(36) 



where the first term represents the homogeneous solution with a prefactor rji which is de- 
termined by the initial conditions at t — 0. The second term represents the inhomogenoeus 
solution which vanishes for t = 0. 

From these solutions of the decoupled problem, one can calculate the corresponding 
solutions of the coupled set of differential equations by x = S ■ y. Finally, imposing the 
initial conditions 



a} a) 



ni, 



t=o 



0, and ((6 f a) - (a f 6)) 



t=o 







(37) 



t=o 



which means that at t — oscillator A has the temperature of its heat bath, oscillator B has 
zero temperature (which is the temperature of the second heat bath), and that at t — the 
photon transfer is just turned on, one can determine rji. The solution for x fulfilling these 
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initial conditions can be written as 



EiWe^ + ^fe^-l] 



(38) 



where 



and 



9 



T)\ = Tlx 

V2 = ni 
V3 = ni 



2k,i + \i 



2/3 — 2/2 



and yi 
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2k,2 + Aj 





N 


yi 


-J/3 




N 


1)2 


-2/i 



N 



N = (x 2 - £i)(yi - y 3 ) + (x 3 - xi)(j/ 2 - yi). 



(39) 

(40) 
(41) 
(42) 
(43) 



In Figs. [3] we show some plots of the analytical expressions for (a^a)/ni, {b^b)/ni, 
Im^tfa) — (a' ( b))/ni and over time choosing K\ = k,2 = K and different coupling strengths 
g/n between oscillator A and oscillator B. First of all, we can see that the mean occupation 
number (a^a) of oscillator A reaches in the long-time limit a constant value smaller than 
its initial value ni, whereas the mean occupation number (tfb) of oscillator B converges in 
the long-time limit to a constant value larger than its initial value 0. Therefore, we ob- 
serve an energy transfer from oscillator A to oscillator B. This energy transfer is related 
lm((tfa) — (a^b))/ni as we will see later. For small coupling strengths the evolution of the 
occupation numbers is monotonic, but for large coupling strengths one can observe an os- 
cillatory behavior. Hence, the energy is going back and forth between both oscillators until 
the steady state is reached. This behaviour is analogous to the Rabi oscillations observed 
for a two-level system coupled to a field including radiative damping, two coupled qubits, or 



coupled exciton-plasmon systems for instance [l|, I22l-l24| . In such systems Rabi oscillations 



are found in the strong-coupling limit which corresponds in our case to g/n 3> 1, whereas 



22 



24]. 



in the weak-coupling limit (here C 1 ) these oscillations are damped out 

In Fig. H]we show Im((6'a) — (a^b)) jri\ over time for g/n = 10 together with (a^a) jri\ and 
(tfb)/ni. As we will see later, Im((tfa) — (a)b))/rii is proportional to the energy transfer 
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Figure 3: Plot of (a) (a^a)/ni, (b) (tfb)/ni and (c) Im((¥a) — {a)b))/n\ over time for k\ = K2 = k 
and different coupling strengths g/n. For large t the solutions converges to the steady state solutions 
found in Eqs. (|23p and (|24p for Q a b = ri{, a = — (^l + ^2) and H2 = 0. 



between oscillator A and B. From Fig.|4]we see that the oscillations of (b'b) are phase shifted 
by 7r with respect to the oscillations of (a^a) and the oscillations of lm((tfa) — (a^6))/ni 
have a phase shift of n/2 with respect to (a^a) and (b^b). 

Since we are interested in the energy transfer from oscillator A to oscillator B let us 
further analyse (b^b) and (b^a) — (a^b) and derive the short- and large-time limits. From the 
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Figure 4: Plot of Im((^o) — (atb))/n\, (a)a)/n\, and (frb)/n\ over time for ki = K2 = 1 = k and 
0/k = 10. 



analytical solutions we find 



2+2 



(6 f a) - (a ] b) 



^1 ' ^2 



1n{\gt 



t < 



(44) 



(45) 



(K1+K2) 1 I 9 2 
K1 K 2 

We see that (tfa) — (atb) is linear in gt and (tfb) is quadratic in gt for times £ <C , . 
For times i ^> k^ 1 , 1 we retrieve the steady-state results derived in Eqs. fl23|) and fl24"|) for 
the case f2 a ft = f2& a = — («i + k 2 ) and yT 2 = 0. Furthermore, we see that for t 3> ftj -1 , k^ 1 

(6+6) oc g 2 (l 



(b ] a) - (a+6) oc 0( 1 



^1^2 
^1^2 



(46) 
(47) 



From this one can expect that (b'b) has only even orders in g, whereas (b*a) — (a'b) has only 
odd orders in g. 



C. Perturbation expansion 

To complete the analysis, we also derive the perturbation result of the set of differential 
equations (jl7p - (l2Up for (b'b) and (b'a) — (a^b) in orders of g. Here again we first consider 
the general case Q ab ^ Q ba . Expanding (at a) = (a*a)W + (ata)^ + (a)a) { - 2 ' ) + (tfb) = 
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+ (tfb)^ + (tfb)( 2 > + ... etc. and solving the dynamical equations for the different 
perturbation orders. We find for the lowest orders giving nonvanishing contributions 



(6+a) (1) - (atfe)« 



10711 



n 



ab 



n 



and 



(2) 



+ 



6a 



-2K 2 t 



(48) 



in 



a6 



6o 



_n a f,(2fi;2 + fi a fc) n6 a (2K2 + fift a ) 2k 2 n^f^ 

As expected (6+6) W = 0. 

To compare these results with the analytical solutions we set u a = Ub = u) sp , i.e. Q a b 
Qba = —(ki + ^2)- Then we find the short and long time limits 



(49) 



(6t a )(D - (at 6 )(D 



(50) 



and 



(2) 



riig 2 t 2 fCKj 1 ,^ 1 



9 "1 



(51) 



(ki+k 2 )k2 

By comparing these perturbation results with the corresponding analytical expressions in 
f !44p and ( |45p we can infer that in general the perturbation results are valid for g 2 <C k±K2- 111 
particular, we find that for t <C k^ 1 , 1 the perturbation result coincides with the analytical 
solution. 



IV. HEAT TRANSFER RATE 

Before we apply our model to describe the dynamics of heat transfer between two nanopar- 
ticles, we will define the energy transfer rate and discuss under which circumstances it can 
be derived from Fermi's golden rule. 

A. Defining the energy transfer rate 

For two classical harmonic oscillators A and B the transferred energy per unit time from 
A to B can be described by the rate of work done on oscillator B by oscillator A 



P = k(x B - x A ) ■ x B 



(52) 
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where k is the spring constant for the spring between the oscillators and xa and xb are the 
displacements of the oscillators from their equilibrium position. The quantum mechanical 
analog can be derived from that expression by replacing 



x A ^(a! + a)J— ^— , (53) 



2mw sp 



^-►(ftt + ft) / (54) 



2mw sp 



x A = — -»> i(a f ~ a ) a/tt^' ( 55 ) 
m V 2m 



± B = ?*-+iV- b )J%*. (56) 
m V 2m 

Keeping only terms for which a photon is interchanged between both oscillators, i.e. which 
allow for the transition processes |ni,n 2 ) — > \fii — l,n 2 + 1) and |ni,n 2 ) —> \n\ + l,n 2 — 1), 
we obtain 

P — ki (atf - ba)). (57) 

In order to relate the prefactor to the coupling constant g or I used in our model, we start 
with the classical expression for the potential energy 

H^^xa-xb) 2 (58) 

and replace the displacements by the expressions in Eqs. ( I5"3"j) and (]54p . Keeping again only 
terms allowing for photon exchange between the oscillators we find 

Hj = -k — - — (a)b + b ] a). (59) 
2muj sp 

By comparing this result with the interaction Hamiltonian of our model in Eq. (J5]) we find 
that g = —k/2mu) sp . Hence for the transferred power 

P = hu sp (-ig)(atf - ba ] ). (60) 



1. Full expression 

Since we know the analytical solution for (atf) — (a^b) we can immediately write down 
the analytical solution for the mean rate of energy transfer 

(P) = ^ sp (-i<?) [(flfet) - (atft)] = hw sp R (61) 
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where the energy transfer rate R is by means of Eq. (I38p given by 

fl=-i0£je** + ^(e*'-lY|. (62) 

j L * \ /J 

From this expression we can derive the short-time and large-time limits 

{2n 1 gH t^Ki 1 ,^ 1 
2g2 _ i _! _r ( 63 ) 

(*?+T 2 ) 1+ V f»«i ,k 2 
K 1 K 2 

Note that in the short-time limit the rate R is given by d{b^b)/dt taking (Wb) from Eq. (1441) . 
On the other hand, in the long-time limit this observation does not remain true. In the 
regime of Rabi oscillations (g ^> K1/2) the rate of transfer is not a meaningful quantity and 
we have to work with the integrated power. 

In the general case (rii > and ri2 > 0) the expression for the transferred power in 
Eq. ( 152]) has to be augmented by a term — k(xA — %b) • &A quantifying the rate of power 
done on oscillator A by oscillator B. The resulting expression for P and R due to that term 
is just the same as in Eq. ( 162|) . but with n\ replaced by — n 2 . Hence, the overall transferred 
energy has a prefactor (rii — n 2 ) determining the direction of the energy transfer or heat 
flux, uniqely. It follows that in steady state the energy is always transferred from the hotter 
to the colder nanosystem. 

2. Perturbation result 

The corresponding expression for the heat transfer rate for small coupling strength is 
obtained with the help of Eq. f!5Up . giving 

R*-ig[{tfa)W-(cfb)M] « I . (64) 

[ (M+I2) *>«1 1 >«2 1 

As noted above for times t -C aci, «2 the perturbation result and the analytical result coincide, 
whereas for times t 3> Ki, &2 the perturbation result is only valid for g 2 <C k±K2- 

B. Fermi's golden rule 

Finally, we derive the heat transfer rate by using Fermi's golden rule 

Rfgr = ^\{f\Hi\i)\ 2 8(uf - Ui). (65) 
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In our case the initial state is given by a) \i) = \rii,n 2 ) and |/) = \n\ — l,n 2 + 1) describing 
the photon transfer from oscillator A to B and b) \i) = |ni,ri2) and |/) = |ni + l,n 2 — 1) 
describing the photon transfer from oscillator B to A. Inserting the interaction Hamiltonian 
H\ from Eq. (J5]) and taking the difference of process a) and b) we find 

#fgr = 2wg 2 [n 1 (n 2 + 1) - n 2 (n 1 + l)]5(u f - tOi) 

(66) 

= 2ixg 2 [ni - n 2 \S(u}f - 
For the special case that rii = ni, n 2 = n 2 = 0, ojf = Uf,, and Ui = u a we obtain 

R FGR = 2iig 2 n 1 5(uj a - u b ). (67) 

Now, let us take into account that the oscillators A and B have a finite linewidth k\ and 
k 2 through their coupling to the reservoirs. Assuming a Lorentzian profile for the linewidth 
and averaging over these profiles, we obtain (setting uj a = cui, = u sp ) 



Rfgr = 9 nxlix I dui du 2 - -— — k- -z- — ^6{oj 1 -u 2 ) 



{u sp - Wi) 2 + k\ {u sp - u 2 ) 2 + k\ 



2g 2 n 1 



(68) 



Ki + K 2 

Comparing this expression with the perturbation result for the heat transfer rate in Eq. f !64p 
we see that it is valid for t 3> aci, k 2 and g 2 k\k 2 . 

To have a better understanding of the Fermi golden rule expression for the heat transfer 
rate, let us see how the golden rule is derived from (tfb). First, the golden rule is based on 
time-dependent first order perturbation theory in the transition amplitudes. Therefore, for 
(b'b) we have to make a second-order perturbation expansion with respect to the coupling 
strength. The resulting expression for (Wb) is already given in Eq. (H9|) . If we cut off the 
heat baths by setting k\ = k 2 = 0, then we arrive at 

sin 2 ^ A 



(bW 2) =9 2 ni ,\V . (69) 



> 2 ' 

A 1 



where A = uo a — w&. From such expression Fermi's golden rule is usually derived by consid- 



ering the long-time limit 



22] 



R ¥GK =lim^ = g^riiA) (70) 
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which is Fermi's golden rule; we have used that 

sin 2 ff 

lim = 27rt5(A). (71) 

A 

2 

Hence, the heat transfer rate can be derived from (b'b) or Fermi's golden rule when first 
cutting off the heat baths first. The resulting rate is only valid in the long-time limit for 
t ^> A -1 and for small coupling strength, since it is based on perturbation theory. The effect 
of the coupling to the heat baths is taken into account by averaging the rate -Rfgr over the 
linewidth profiles. For two nanoparticles with u a = Ub = oo sp having a finite linewidth K\ 
and K2 the long-time limit condition then translates into t S> 

V. HEAT EXCHANGE BETWEEN TWO NANOPARTICLES - 
DETERMINATION OF THE PARAMETER g IN THE MICROSCOPIC 
INTERACTION HAMILTONIAN Hi 

Up to here we have discussed our model itself in some detail. To relate it to the heat 
transfer dynamics for the resonant radiative heat exchange between two nanoparticles, we 
still need to determine the coupling g for describing the heat transfer rate. To do this, we 
consider the heat flux between two nano-particles in the steady-state and determine g from 
that by comparing the heat transfer rate with the steady-state solution of our model. 



A. Heat transfer rate in steady state 



The heat transfer rate between two nano-particles in steady state was derived for example 



in Refs 



11 



the second kind 



17]. Its derivation is usually based on the fluctuation-dissipation theorem of 



HE 



5[ which serves as the basis of Rytov's fluctuational electrodynamics 



Within this framework the power P transferred between two identical nanoparticles with 



temperatures T\ > and T 2 = is 



111 







An 3 



huriilm(a) — 



(72) 



This expression is only valid for particles with a radius r smaller than the dominant thermal 
wavelength and an interparticle distance d larger than their radii. Furthermore, multiple- 
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interactions are neglected as well as the contribution due to eddy currents. 

Now, since we are interested in the power transferred by the resonant interaction of the 
localized surface modes of the nanoparticles we expect that the main contribution comes 
from the resonance of the polarizability a at frequency u sp which allows for approximating 
the above expression by 



00 dw. 



P « hu^FM I ^lm(a) 2 (73) 







introducing the function 



F(u) 



Therefore the heat transfer rate in steady state is 

-R s t-st = niF(o; S p) J -^lm(a) 2 
= n 1 F(u sp )(Anrl) 



(74) 



2 1 f + °° du f3lm{e) 



(75) 



2 7_oo 4tt3 

To further proceed we assume that we have metallic nanoparticles which can be described 
by the Drude model in Eq. 03]). Inserting this expression in the integrand of i? s t-st we obtain 

,^lr 2 / 3 V r + °° . 1 



u b r 2 { 3 \ f +OQ 



u 2 - oo 2 p + 2iooT) 2 (uj 2 - co 2 p - 2icoT) 2 

(76) 



introducing the surface mode frequency oo sp = lo p / V e <x> + 2 and its damping T = 7/2. 
The frequency integrand has four second-order poles in the complex frequency plane. We 
compute the integral by taking the two poles in the upper half plane into account and obtain 
assuming that u sp 3> T 

oo 2 r 6 / 3 \ 2 

fl „^ lfimp) (_L_). (77) 

B. Comparison with master equation approach 

Since we have the steady-state result for the heat transfer rate between two identical 
nanoparticles in Eq. (177|) we can determine the coupling constant g by comparing it to 
the steady-state solution of our model. Since expression fT72|) does not include multiple 
interactions between the nanoparticles it is in fact a first order perturbation result in a 2 . 
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Therefore, we compare it with the perturbation expression for the steady-state solution in 
Eq. (JMD setting Kl = k 2 = T 

R = |W (78) 

By comparison to i? s t-st we can identify 




Since to leading order for small distances F(u sp ) oc 1/d 6 we have g oc 1/d 3 which is typical 
for a dipole-dipole interaction. 



C. Heat transfer dynamics for two nanoparticles 

Finally, we have all the ingredients to determine the dynamics of the heat transfer rate 
for two nanoparticles supporting surface modes. When assuming that we have two gold 



nanoparticles then ui p = 1.4 x 10 rad/s, = 3.7 and 7 = 2.79 x 



O^s" 1 (see Fig. 2.1 



on p. 25 of Ref. [27] for a fit to the Christy- Johnson data from Ref. [26]). It follows that 
the surface mode resonance frequency is u sp = 5.86 x 10 15 rad/s with the linewidth T = 
1.4 x 10 13 s _1 . If we consider a nanoparticle radius of r$ = 20 nm and an interparticle distance 
of d = 800 nm, and d = 1200 nm we find the heat transfer rate R/T plotted in Fig. [5j For 
times larger than 1/T the heat transfer rate coincides with its steady-state expression. On 
the other hand, in the transient regime one can observe the dynamics of the heat transfer 
rate for times smaller than 1/T « 72fs. As can be seen, for the considered distances we are 
mainly in the weak-coupling regime for which we have derived the coupling constant from 
the known steady-state heat flux expression. For distances smaller than d = 800 nm the 
coupling becomes stronger. Finally, in the strong coupling regime the energy transfer rate 
is not a meaningful quantity so that we consider for this regime (tfb) jH\ rather than R/T. 
In Fig. Owe plot (tfb)/ni for d = 100 nm which corresponds to g/T = 5.8 showing clearly 
the Rabi oscillations. 



VI. CONCLUSION 

We have introduced a model of two coupled dipoles which are both coupled to their own 
heat baths of independent oscillators. This model allows us to describe the heat transfer 
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0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 

tr tr 

(a) d = 800 nm, g/T = 0.6 (b) d = 1200 nm, g/T = 0.4 

Figure 5: Heat transfer rate R/T from Eq. (|62p between two nanoparticles normalized to the 
linewidth T of the nanoparticles surface mode resonances over time in units of T -1 . We use the 
coupling constant g from Eq. (|79|) with tq = 20 nm and d = 800 nm and d = 1200 nm. The 
horizontal line is the steady-state expression from Eq. (|63p . 
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tr 

Figure 6: Mean occupation number {b'b)/ni of the surface mode of nanoparticle B from Eq. (|38p 
over time in units of T" 1 . We use the coupling constant g from Eq. (|79p with r$ = 20 nm and 
d = 100 nm, hence g/T = 5.8. The horizontal line is the steady-state expression from Eq. (|44D . 

dynamics between two nano systems. In particular, we have used this model to determine 
the heat transfer rate between two spherical nanoparticles which is due to the resonant 
interaction of the localized surface modes. Within our model we have derived the expression 
in the limit of short times as well as in the steady state regime and made a connection to 
the usual heat transfer calculations which are based on the fluctuation-dissipation theorem 
considering the steady-state regime only. We have shown that for small coupling strength 
between the nanoparticles the steady-state heat transfer rate can also be determined by 
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Fermi's golden rule when taking into account the finite widths of the resonant surface modes. 
Finally, we have studied the dynamical heat transfer rate for two gold nanoparticles in the 
transient regime for times shorter than the relaxation time of the surface modes and predict 
Rabi oscillations for the mean occupation number of the surface plasmons 
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